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We present a new method for obtaining matched solutions of the rms envelope equations. In this 
approach, the envelope equations are first expressed in Hamiltonian form. The Hamiltonian defines 
a nonlinear mapping, M, and for periodic transport systems the fixed points of the one-period map 
are the matched envelopes. Expanding the Hamiltonian around a fiducial trajectory, one obtains a 
linear map, M, that describes trajectories (rms envelopes) near the fiducial trajectory. Using M 
and M we construct a contraction mapping that can be used to obtain the matched envelopes. The 
algorithm is quadratically convergent. Using the zero-current matched parameters as starting values, 
the contraction mapping typically converges in a few to several iterations. Since our approach uses 
numerical integration to obtain all the mappings, it includes the effects of nonidealized, z-dependent 
transverse and longitudinal focusing fields. We present numerical examples including finding a 
matched beam in a quadrupole channel with rf bunchers. 



INTRODUCTION 



In 1959 Kapchinskij and Vladimirskij published the first envelope equations governing two dimensional beams with 
space charge [1]. Although they assumed an unusual phase space distribution (a 6-function whose argument was a 
linear function of the Courant-Snyder invariants) their result was important because it described, in terms of a set of 
ordinary differential equations, the self-consistent transport of finite emittance beams in strong focusing systems. More 
than 10 years later Sacherer and Lapostolle separately showed that, for beams with elliptical symmetry, one could 
derive rms envelope equations that were satisfied by beams in general and not just KV beams [2,3]. They showed 
that for linear external fields (but no such restriction on the beam self-fields) one could obtain a set of equations 
that involved only second moments, but to achieve this they had to allow the beam emittance to appear in the rms 
equations as an unknown function of time. Sacherer also showed that one could derive envelope equations for three 
dimensional beams (i.e. a six dimensional phase space) assuming ellipsoidal symmetry; unlike the two dimensional 
case, the equations depend on the precise form of the distribution through a parameter called A3, but that dependence 
turns out to be very weak. Through the years activity has continued in the area of envelope equations. Motivated by a 
desire to analyze the transport of intense electron beams in gas, Lee and Cooper derived rms equations for cylindrically 
symmetric beams including scattering; additionally, their formulation included nonzero canonical angular momentum 
and acceleration [4]. In order to model beam transport systems and rf linacs, Crandall developed the now widely 
used codes TRACE and TRACE3D [5]. More recently, Struckmeier derived envelope equations starting from the 
Fokker-Planck equation [6]. He used the formalism to estimate emittance growth due to intra-beam scattering in 
storage rings. 

In this paper we will emphasize the application of envelope equations to quadrupole channels and rf linacs. Envelope 
equations have turned out to be extremely useful in the early stages of the design of these systems. Our approach 
allows a very accurate treatment of the beamline elements, since it takes into account ^-dependent longitudinal and 
transverse fields of rf gaps and focusing magnets. The fields can be in the form of measured field data or analytic 
functions that approximate field data. 

Often one is interested not only in numerically integrating the rms equations but also in finding initial conditions that 
result in periodic solutions when the beam transport system is itself periodic. In this paper we will present an efficient 
method for finding these matched solutions. Our method rests on the fact that the rms equations can be represented 
as a Hamiltonian system. While this system's nonlinear behavior can be found through numerical integration of the 
rms equations, its linear behavior can be obtained using standard techniques in accelerator physics for computing 
linear maps. Together, the linear and nonlinear maps can be used to construct a quadratically convergent iterative 
procedure for finding matched rms envelopes. Typically, we have found the procedure to converge in a few to several 
iterations, even under conditions of extreme space charge depression. 



I. TWO DIMENSIONAL SYSTEMS 



A. Overview 



Consider a particle beam propagating in a quadrupole channel. Suppose that the beam is long compared with its 
transverse dimensions, and that we can neglect any longitudinal variation when calculating the beam self-fields. We 
will neglect image charge effects, and we will suppose that the beam is launched along the axis of a perfectly aligned 
transport system. We will use the longitudinal coordinate, z, as the independent variable. The canonical coordinates 
and momenta for the transverse phase space are denoted (x,p x ,y,p y ). Let the vector potential associated with the 
quadrupoles be given by 

A x =A y = 0, (1) 
A z = l -g{z)tf - x% (2) 

where g(z) denotes the magnetic quadrupole gradient. Let ^ self denote the scalar potential associated with the 
self-fields and, neglecting transverse currents, suppose that the associated vector potential is given by 

A x =A y = 0, (3) 

Az = ^* self , (4) 
c 

where (3 c is the velocity on the design trajectory. Rather than working with the variables (x,p x , y,p y ) it is convenient 
to define dimensionless variables (x,p x ,y,p y ) according to 

x = x/l, p x = p x /p , (5) 
ti = y/l, Py=Py/Po, (6) 

where p denotes the momentum on the design trajectory (i.e. p = y (3 mc) and where I is a scale length [7]. The 
Hamiltonian (in MKSA units) governing these variables is given approximately by 

H{x,p x ,y,p y ,z) = ^{Pl +P 2 y ) 

+ l -^(x 2 -?) + ^(lx,ly,z), (7) 



where 



and where K is the generalized perveance, 



Also, i& is related to ij self according to 



K z ) = {llPo)g{z) 
2ire p v z Y 



fself = (10) 
47re V 1 

where A is the charge per unit length measured in the lab frame, A = I/v . Note that we have expanded the relativistic 
Hamiltonian to second order in the phase space variables, with the exception of the scalar potential, as is the standard 
procedure for deriving rms envelope equations. For the remainder of the discussion of two dimensional systems we 
will set !=lm. 

Following the usual procedure one can obtain equations for the rms envelopes, X and Y [2,3]. The envelope 
equations are given by 

d?X , v K/2 El n 

^ + ^-yT7-# = ' 

^-*y--^-4 = o, (ii) 

dz 2 X + Y Y 3 K ' 
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where £ x and £ y denote unnormalized rms emittances. Since these are rms equations the factor in the space charge 
term is K/2, whereas it would be 2K for the KV equations. The envelope equations are derivable from a Hamiltonian, 
H env , where 

H«"{X, P X ,Y, P y ) = l -{Pl + ||) + l -{P* y + ^) + k -{X* - Y 2 ) - (if/2) log(X + Y). (12) 

Lastly, one can use the envelope Hamiltonian to define depressed phase advances, \i x and (j, y , of a particle in an 
equivalent KV beam. Referring to the envelope Hamiltonian, we regard the phase advances as coordinates and the 
emittances as a canonically conjugate momenta [8]. (Since the phase advance appears nowhere in the Hamiltonian, 
the emittance is constant, as expected). Taking this view, we obtain 

-F f 

/dz 

For a KV distribution, these formulas apply to all the particles in the beam (since the forces are linear). For other 
distributions they apply to the equivalent KV beam (i.e. a KV beam with the same rms values). The phase advance 
equations can be integrated along with the envelope equations. 



B. Constructing the contraction map 



Consider a periodic transport system with period L. Let £ = (X, P X ,Y, P y ). As shown above the envelope equations 
are derivable from a Hamiltonian; hence they define a symplectic nonlinear mapping M that maps initial state vectors 
into final state vectors: 

£ fin = MC n . (14) 

If we consider transport through one period of the transport system, then a matched envelope is simply a fixed point 
ofM: 

M( = (. (15) 

Techniques for finding fixed points of symplectic maps are widely used in accelerator physics [9]. For example, they 
are used to find the off-momentum closed orbits of particles in circular machines. In our case, however, we will use 
the techniques to find the fixed points of an envelope map, not a particle map. The approach is based on the fact that 
the machinery exists to compute the action of the nonlinear map M as well as its linear part, M. This makes it easy 
to construct a contraction map based on a Newton search procedure to find the fixed point. As a result, the method 
is quadratically convergent. As an illustration, consider the problem of finding solutions of the equation g(x) = x, 
or equivalently, of finding roots of the function f(x) = g(x) — x. The Newton search algorithm defines a contraction 
mapping C that, for sufficiently close starting values, converges to a root [10]. If x n is the value of x on the nth 
iteration, then applying the contraction map C to x n produces a value at the next iteration: 

= Cx n = x n - 4^r = * n ~ ~ • (16) 

/'(a") l-g'{x n ) y 1 

It is easily shown that if x n is within e of a root then Cx n deviates by an amount proportional to e 2 . 
For a multidimensional system, the contraction map is given by [9] 

C +1 = CC =C -{I- M) -1 (£ n - MC), (17) 

where J is the identity matrix and M is the matrix associated with linear part of M. To obtain M, we first need 
to linearize H env about a "given" fiducial trajectory, £ g . Let £ = C ~ Cg- The quadratic part of the Hamiltonian 
governing these deviation variables, which we will denote H2, is given by 

1,-2 "2, X\, 3£ 2 (K/2) , Y\ , 3£ 2 (K/2) , (K/2)XY , x 



3 



The equation of motion for M is well known [11]: 



^ = JSM, (19) 
dz 

where the symmetric matrix S is related to H2 by 

1 4 

a, 6 = 1 

and where the matrix J is the fundamental symplectic two-form (i.e. the nonzero elements of J are given by 
J12 = J34 = 1, J21 = J43 = — 1)- Comparing Eq. (18) and Eq. (20) (with £ = (X , P x , Y, P y )) one can immediately 
identify the matrix elements of S. That is, Su is the coefficient of \X 2 , S22 is the coefficient of ^P 2 , Sis = £31 is 
the coefficient of XY , etc. 

Summarizing, Eq. (17) defines a contraction map for finding matched rms envelopes. In order to evaluate the right 
hand side of the equation, one must use numerical integration to compute the following: (1)A4£, which is just the 
numerical solution of Hamilton's equations with the Hamiltonian of Eq. (12); and (2)the matrix M , which is obtained 
by numerically integrating Eq. (19). These quantities are computed at every iteration until the difference between CC, 
and £ is sufficiently small. In the calculations below we consider the map to have converged when |£ — C£|/|£| < 10 -8 . 
We use the zero current matched values as starting values. For systems where these cannot be found analytically, we 
begin by integrating the envelope equations with zero current and zero emittance; this is equivalent to computing the 
beta functions, from which we obtain the zero current matched envelopes. 



C. Example: matching in a quadrupole channel 



As a numerical example, consider a beam in a magnetic quadrupole channel of the "FODO" type. The quadrupoles 
in this example are idealized as having a constant gradient over their length, though since numerical integration 
is used anyway, it would be little or no extra work to use analytic models or numerical values for the quadrupole 
gradients. The period length L = 24 cm and each quadrupole is 6 cm long. The channel is designed to transport a 
beam of 10 MeV protons having rms emittances £ x = £ y = 1 x 10 -6 m-rad. The zero current phase shift per cell is 60 
degrees, which requires g = 78 T/m. At zero current, the matched rms beam sizes are given by X = 0.6310 mm and 
Y = 0.3796 mm. With these as starting values the contraction map converges in 5 iterations when the beam current 
is / = 8.4 amp, as is shown in the following computer output listing: 
Enter current, x-emittance, y-emittance: 
8.4 l.e-6 l.e-6 

Zero current matched rms values: 
X=6.3103d-04, P_x=-2.0322d-19 
Y=3.7961d-04, P_y=-3 . 3781d-19 
Zero current phase shifts per cell: 
sigma0_x=59.971, sigma0_y=59 . 971 



Starting contraction mapping. . . 



iteration 1 
iteration 2 
iteration 3 
iteration 4 
iteration 5 



delta= 
delta= 
delta= 
delta= 
delta= 



3.326906d-00 
8.062686d-01 
3.867283d-03 
1.059674d-04 
7.206193d-10 



search converged 



Matched rms values: 
X=1.9492d-03, P_x=l . 45795d-10 
Y=1.2217d-03, P_y=2.42433d-10 
Phase shifts per cell: 
sigma_x=6 . Oil , sigma_y=6 . Oil 



Note that the convergence is very rapid even though the depressed phase advance is only 6 degrees per cell. 
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II. THREE DIMENSIONAL SYSTEMS 



A. Overview 



Consider a particle beam propagating in a beamline consisting of quadrupole magnets and cylindrically symmetric 
rf cavities. These elements make it possible to approximately model structures such as drift tube linacs and coupled 
cavity linacs. The potentials associated with the quadrupole magnets and the beam self-fields are the same is in the 
two dimensional case (but we will use the notation g m instead of g to denote the magnetic quadrupole gradient). 
Additionally, the vector potential associated with an rf cavity is of the form 

e'(z) 

A x = -^-xsm(w a t + 6) 
2,ui a 

e'(z) 

A y = -^ysm(u a t + 0) 

A z = ~—{e{z) - ^[e"{z) + ^fe{z)]}am{u, a t + 6), 
UJ a 4 C 

(21) 

where the electric field at r = is given by 

E z {r = 0) = e(z)cos(u a t + 6), (22) 

and where a prime denotes d/dz. 

As before, we use the longitudinal coordinate z as the independent variable. Now there are six canonical coordinates 
and momenta, denoted (x,p x ,y,p y ,t,p t ), where t denotes a particle's arrival time at the location z and where its 
canonically conjugate momentum p t is the negative of its total energy. As above, we will define variables that are 
dimensionless deviations from the design trajectory. First, let (t ,pt ) denote the design trajectory (along with 
x = p x = y = p y = 0). Using the above potentials, it follows that the equations of motion for the design trajectory 
are given by 

, = -Pto/C 

o 1—2 FT v ' 

VPto - m C 

p' t0 = -qe(z)cos(u a t + 6). (24) 

Note that, in the two dimensional case, the momentum on the design orbit was a constant, and it was natural to scale 
the transverse momenta by this quantity. But in the three dimensional case, the rf fields can accelerate the beam. 
Thus, we will scale the transverse momenta by a parameter 8 which is unrelated to the design momentum. (In fact, 
later it will be convenient to set 8 = mc). The transverse coordinates are scaled by a parameter I. Lastly, the time t 
will be scaled by a quantity ui, so that the times are really phases. Often one would choose ui = ui a , but due to the 
frequency changes typical of proton or ion linacs, there are instances where one would choose nui = ui a , where n is an 
integer. In summary, the dimensionless deviation variables are given by 

x = x/l, p x = p x /8 (25) 

y = y/i, Py =Py/S (26) 
t = w(t-t ), p t = {p t - Pt0 )/{wl8) (27) 

The single particle Hamiltonian, H(x,p x ,y,p y ,i,pt, z), paraxial in the external fields, is given by 

qe' sin d>. , „ m 2 u> 2 l8 , w„oe sin d>. ^, p , - , . ~ . -, . , 

_q_ ^ ( __^ + + _^ p 2 _ a« V. p + £ * Jjj, ly t tu + to , Z ), (28) 

where the synchronous phase, <f> s , is given by 

d> s =u n Uz) + 6. (29) 
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In the above Hamiltonian, ^ is obtained from Eq. (10) with A = I/v , and K is the perveance obtained from Eq. (9) 
using the average current J. Note that if the scaling parameter ui is chosen to be the frequency of the bunches, then 
the charge per bunch Q is related to J according to Q = (2-k/ui)I. In what follows, it will be convenient to define a 
quantity u , where 



(30) 



X" 


+ P± X ' 




19rS x 




Po 


Po 


Po 


Y" ■ 




19m y 


99rS Y 




Po 


Po 


Po 


T" - 


f 3— T' 

Po 


qu a /c 2 
llPoPo 


e sin <f> s T 



where 



7oPo 

u = — rr . 

Following Sacherer [2], one can obtain equations for the rms envelopes, X, Y and T: 

-^XG 31l{ X,Y,u o T)-(±f-^ = , 

-^TGMX^UoT)-^)^^, (3!) 

1 u) a (w/w a ) 2 , 

9rf = -[^-esin^ s e cos<£ s J. (32) 

z c z v 

In the above equations, £ njX , E n ^ y and £ njt denote normalized rms emittances. The quantity A3 is a geometrical factor 
that depends on the details of the charge distribution within the bunch, but as Sacherer pointed out it is not very 
sensitive to the details and has a value approximately equal to l/(5\/5) for a wide variety of distributions. Lastly, 
the quantity G is a space charge term defined by 

- 3 f°° ds 

G mnp (x,y,z)--J^ (a , 2+s)m /2( y 2 +s) n/2 (2 2 + s) p/2- ( 33 ) 

Note that these rms equations are not expected to accurately model the bunching process. The reason for this is 
twofold: (l)By our paraxial expansion of the single particle Hamiltonian, we assume that the external rf fields vary 
linearly across a bunch; (2)The space charge terms are based on the fields of an isolated bunch of charge, not a train 
of bunches. 

Lastly, the rms equations are derivable from the following envelope Hamiltonian: 

H- { Z, P., Y, P„ T, P.) = J-(P? + |f ) + JL(p ; + |f ) + ^(P? + |f) 
,1 l 9m (Y 2 Y 2^ 9kr S , 2 , y 3N g^ e sin (j> s 2 Ku 2 ttX 3 

+ ^p7 {x ~ y) ~ ~^t {x +y) ~ 2^16 T + — r~ Glll(x ' y ' u ° n (34) 



B. Constructing the contraction map 

As in the two dimensional case, we need to compute the matrix M that describes the linear behavior of the system 
governed by H env . This in turn requires that we know the matrix S, which appears in Eqs. (19) and (20). Linearizing 
H env around a given fiducial trajectory £ g we obtain the following nonzero matrix elements of S: 



6_ 

lp ' lp u 2 " 



322 — ^44 — — , See — j, (35) 

_ 1 



S\\ = -J\9m - 9rf) + — ^ 1 {3X g G 511 - G311), (36) 

S 33 = -j\-9m - 9rf) + ]^~^T~ ^ 1 y 3Y g ° 151 ~ G ^i), (37) 
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'55 



qu a esm4> s 8 3£ 2 t if ^7^3 



3^T 2 G 115 - G 113 ), (38) 



5l3 — 531 — ^ -^5^5^331, (39) 

5i5 = 55i = ^ XgTgGsis, (40) 



'35 — ^53 



— ^T,G 133 . (41) 



Note that we have used the notation G mnp to denote G mnp (X, Y, u T) in the above equations. 

The following section contains a numerical example of finding matched beams in a three dimensional system. But 
before continuing we need to specify our choice of scaling parameters: 

nui = u) a , (42) 
6 = mc, (43) 
I = c/u>. (44) 

In the above equations, ui is simply the frequency of the bunches, normally a harmonic of the rf frequency. If every 
rf bucket is filled, n = 1; if every other bucket is filled, n = 2; and so on. By choosing uil/c = 1, it means that X and 
Y need to be multiplied by the inverse wavenumber, & _1 = c/u>, to convert them to dimensional quantities. Lastly, 
note that with this choice of parameters the phase advances are given by 



(J-x 

fJ, y 
Mi 



&n , x f 

~rj 7o p x 2 ' 

&n , y f 

£ n t f dz 

~t]^w^- (45) 



C. Example: matching in a quadrupole channel with bunchers 

We will consider the same FODO channel as described previously, but in addition rf gaps will be inserted between 
each quadrupole. This is shown schematically in Fig. 1. Though we could use numerical field data or analytic 
functions that approximate field data, for the sake of illustration we have assumed that e(z) is a sum of two identical, 
longitudinally separated Gaussians: 

e(z) = E max (e^- Z ^' 2 " 2 + e^-^l 2 " 2 ), (46) 

with z\ = 6 cm, Z2 = 18 cm, a = 5 mm and E max = 20 MV/m. Each cavity is assumed to have the same frequency 
and phase (see Eqs. (21)): /„ = w a /27r = 361.75 MHz, 6 = 90.06 degrees. With these choices of frequency, phase, and 
cavity separation, a synchronous particle crosses the gaps with a phase of roughly —90 degrees, the net result being 
no acceleration or deceleration. Fig. 2 shows the synchronous particle's energy as a function of z. 

The transverse emittances are chosen to be the same as in the two dimensional example. Since we now need 
dimensionless normalized values, we have to take the numbers from the two dimensional example and multiply them 
by a factor (p/ (3 /l). This factor equals 1.11, and so the resulting dimensionless normalized rms emittances are given 
by £ n ,x = £n,y = 1-11 x 10 -6 . For simplicity this is also the value we choose for £ n ,t- As before, the channel is 
designed to transport a beam of 10 MeV protons. At zero current, the matched rms beam sizes are X = 4.968 x 10 -3 , 

Y = 3.009 x 10 -3 and T = 3.416 x 10 -2 rad. (In physical units, the transverse values are X = 0.6553 mm and 

Y = 0.3969 mm, since I = 0.1319 m). At zero current, the rf defocusing effect of the gaps depresses the transverse 
phase advances from 60 degrees to around 55 degrees per period; the temporal (i.e. longitudinal) phase advance is 
slightly more than 30 degrees per period. At a beam current of 150 ma, the contraction map converges in 6 iterations: 
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f req, curr ,x-emit ,y-emit ,t-emit (normalized) : 
361.75e6 0.150 l.lle-6 l.lle-6 l.lle-6 



Zero current matched rms values: 
X=4.9679d-03, P_x=-2 . 0090d-10 
Y=3.0094d-03, P_y=6 . 1963d-12 
T=3.4158d-02, P_t=-2 . 7831d-12 
Zero current phase shifts per cell: 
sig0_x=55 . 377 , sig0_y=55 . 376 , sig0_t=31 . 434 



Starting contraction mapping. . . 



iteration 


1 


delta= 


9 


460663d- 


01 


iteration 


2 


delta= 


4 


247934d- 


01 


iteration 


3 


delta= 


3 


944055d- 


02 


iteration 


4 


delta= 


5 


912194d- 


04 


iteration 


5 


delta= 


1 


848957d- 


06 


iteration 


6 


delta= 


4 


322956d- 


09 



search converged 

Matched rms values: 

X=8.5011d-03, P_x=1.1776d-09 

Y=5.3076d-03, P_y=2 . 4236d-09 

T=9.9245d-02, P_t=l . 9463d-12 

Phase shifts per cell: 

sig_x=18 . 301 , sig_y=18 . 300 , sig_t=3 . 724 

Fig. 3 shows the matched rms envelopes in the transverse directions, and Fig. 4 shows the matched temporal envelope. 



SUMMARY 

The purpose of this paper has been to present a method for finding matched rms envelopes in beamlines consisting of 
quadrupole magnets and cylindrically symmetric rf gaps. We did this by casting the envelope equations in Hamiltonian 
form and constructing a contraction map based on a Newton search procedure. Our numerical examples show that 
the map converges very rapidly even under conditions of significant space charge depression. 
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FIG. 1. One cell of a quadrupole channel with rf bunchers, showing quadrupole gradient (solid line) and longitudinal electric 
field (dashed line). 
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FIG. 2. Energy versus z for a synchronous particle propagating in one cell of the quadrupole channel with bunchers. 



9 



0.008 



co . 0075 



0.007 



X 0.0065 



0.0055 



.005 




0.12 0.15 
z (m) 



FIG. 3. Matched transverse envelopes for the quadrupole channel with bunchers. Since X and Y are dimensionless, they 
must be multiplied by the scale length I = 0.1319 m to convert them to meters. 
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FIG. 4. Matched temporal (i.e. longitudinal) envelope T for the quadrupole channel with bunchers. 
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